Para nuestro trabajo hemos elegido los datos de una base de datos pública que hemos encontrado en el sitio web kaggle.com:
https://www.kaggle.com/uciml/biomechanical-features-of-orthopedic-patients (datos recodidos inicialmente de Lichman, M. (2013). UCI Machine Learning Repository [http://archive.ics.uci.edu/ml]. Irvine, CA: University of California, School of Information and Computer Science)
Hemos eligido este conjunto de datos (data-set) porque: - Uno de los problemas frecuentes dentro de la medicina (en el ámbito contemporáneo) es la optimización de recursos donde sea posible. Y un área de aplicación de los análisis de datos es la etapa del diagnóstico de algunas enfermedades en las etapas tempranas. Por ejemplo, el diagnóstico de la hernia intervertebral exige utilizar prubas caras y que tardan mucho en realizarse, cuando hay alternativas que quizás no son tan fiables, pero puedan ayudar a la hora de selecionar pacientes para la siguiente etapa de estudios. En este caso se puede hablar de comparación entre los estudios con tomografía computarizada vs estudios de resonancia magnética. En este sentido, el hecho de definir unas variables posturales pélvicas determinadas y agrupar los pacientes según ellas podría ayudar a este tipo de diagnósticos, y el análisis de los datos de pacientes con patologías determinadas en comparación con personas sin esas patologías es crucial para poder definir unos valores de diferenciación para dicho diagnóstico. - Es un data-set público y fácil de acceder. - Es un data-set relativamente limpio y ya preparado para realizar análisis de datos sobre él.
Este conjunto de datos tiene 2 archivos en formato “.csv” (column_2C_weka y column_3C_weka), que se pueden unir y valorar conjuntamente. El data-set es una base de datos que contiene informacaión sobre 310 personas, en las que se evaluaron la posición de sus vértebras de la parte lumbar utilizando varios parámetros (véase la tabla 1 y una explicación visual del artículo 1). Dentro de esta muestra había 100 personas que consideramos “sanas” o “normales”, sin patologías vertebrales o pélvicas, y 210 (“anormales”) que consideramos que tenían algun tipo de patología: “hernia” o “spondilolistesis”.
tabla 1:
| Variable | Expicación del variable |
|---|---|
| pelvic_incidence, numérica, en grados | Ángulo entre la línea que pasa del centro de la lámina terminal del S1 al centro de las cabezas femorales, y la línea perpendicular a la lámina terminal de S1 |
| pelvic_tilt, numérica, en grados | Ángulo entre la línea que pasa del centro de la lámina terminal del S1 al centro de las cabezas femorales y el eje vertical |
| lumbar_lordosis_angle, numérica, en grados | Ángulo entre la lámina terminal superior del L1 y S1 |
| sacral_sclope, numerica, en grados | Ángulo entre la línea perpendicular a la lámina terminal de S1 y el eje horizontal |
| pelvic_radius, numérica, en grados | Ángulo entre el eje vertical y una línea entre centro de las cabezas femorales hasta el borde posterior de la lámina terminal de S1 |
| degree_spondylolysthesis, numérica, en mm | Distancia de desplazamiento de la vértebra superior sobre la vértebra inferior (positivo: anterior/ negativo: posterior) |
| class | Es posible utilizar esta variable como variable de tipo character o de tipo factor. |
| class_2c, factor | Clasificación de los sujetos como “Normal” y “Abnormal” |
| class_3c, factor | Clasificación de los sujetos como “Normal”, “Hernia” y “Spondylolisthesis” |
En esta imagen presentaremos la visualización de lo ángulos medidos en la base de datos.
Las bibliografías utilizadas para la comprensión de estas variables se han extraído de las siguientes referencias:
1. Sergides IG, McCombe PF, White G, Mokhtar S, Sears WR. Lumbo-pelvic lordosis and the pelvic radius technique in the assessment of spinal sagittal balance: strengths and caveats. Eur Spine J. 2011;20 Suppl 5(Suppl 5):591-601. doi:10.1007/s00586-011-1926-z
2. https://www.orthobullets.com/spine/2071/lumbar-spine-anatomy para la fecha de 25.12.2021
# Vamos a introducir los datos en RStudio. Para esto vamos a utilizar la librería básica de R
column_2C <- read_csv("archive/column_2C_weka.csv",
col_types = cols(pelvic_incidence = col_number(),
`pelvic_tilt numeric` = col_number(),
lumbar_lordosis_angle = col_number(),
sacral_slope = col_number(), pelvic_radius = col_number(),
degree_spondylolisthesis = col_number()))
column_3C <- read_csv("archive/column_3C_weka.csv",
col_types = cols(pelvic_incidence = col_number(),
pelvic_tilt = col_number(), lumbar_lordosis_angle = col_number(),
sacral_slope = col_number(), pelvic_radius = col_number(),
degree_spondylolisthesis = col_number()))
# Vamos a ajustar los nombres de las variables y confirmaremos que son iguales
colnames(column_2C) <- colnames(column_3C)
colnames(column_2C) %>%
set_names() %>%
map(~(column_2C[,.x]==column_3C[,.x])) %>%
as.data.frame() %>%
summary() #utilizaremos la función "map" de la librería "purr" de tidyverse, aplicando la función a todos los elementos de las columnas correspondientes.
## pelvic_incidence pelvic_tilt lumbar_lordosis_angle sacral_slope
## Mode:logical Mode:logical Mode:logical Mode:logical
## TRUE:310 TRUE:310 TRUE:310 TRUE:310
##
## pelvic_radius degree_spondylolisthesis class
## Mode:logical Mode:logical Mode :logical
## TRUE:310 TRUE:310 FALSE:210
## TRUE :100
# Y cambiamos el nombre de la columna "class" para poder juntar los dos data-frames.
column_2C<-rename(column_2C, "class_2" = "class")
column_3C<-rename(column_3C, "class_3" = "class")
# Creamos un data-frame con información de los dos archivos iniciales.
ortho_df <- merge(column_2C, column_3C, by = c("pelvic_incidence", "pelvic_tilt", "lumbar_lordosis_angle", "sacral_slope", "pelvic_radius", "degree_spondylolisthesis"))
# Convertimos la variable "class" en un factor, con los niveles que van de nivel Normal a nivel con alteraciones más faáiles de encontrar (no se puede evidenciar la hernia intervertebral en la radiografía; pero la espondilolistesis, sí).
ortho_df <- ortho_df %>%
mutate("class_2" = factor(ortho_df$`class_2`,
levels = c("Normal", "Abnormal")))
ortho_df <- ortho_df %>%
mutate("class_3" = factor(ortho_df$`class_3`,
levels = c("Normal", "Hernia", "Spondylolisthesis")))
# En algunas ocasiones, para el objetivo de visualizar los datos con mayor facilidad, vamos a crear un data-set con:
ortho_df2 <- ortho_df %>%
pivot_longer(cols = c(pelvic_incidence, pelvic_tilt, pelvic_radius, sacral_slope, lumbar_lordosis_angle),
names_to = "type_of_measure", values_to = "angles") %>%
select(type_of_measure, angles, everything())
ortho_df3 <- ortho_df %>% filter(degree_spondylolisthesis<400) # quitamos del data.frame un valor de la muestra evidentemente erróneo (paciente con un degree_spondilolistesis de 418 mm)
De todas las cuestiones que podemos realizar con estos datos, seleccionaremos las dos que, para nosotros, son las más destacables y extrapolables a una utilidad médica real.
La primera cuestión sería relativa a la comparación de los grupos de estudio. Pregunta 1: ¿Existen diferencias significativas entre las variables angulares del grupo “Normal” respecto al grupo “Abnormal” (personas con hernia o con espondilolistesis)?. En otras palabras, ¿el hecho de tener hernia o espondilolistesis afecta significativamente en alguna de las variables de estudio? Su respuesta nos permitirá discernir qué métodos de medida difieren estadísticamente respecto al grupo “Normal” y llegar a un mejor acercamiento a nivel de diseño de prótesis o clasificación de características de pacientes en bases de datos.
Para resolver esta pregunta compararemos los datos de los grupos a lo largo de este estudio, visualizaremos los estadísticos básicos de nuestra muestra y realizaremos un análisis de varianza (ANOVA) de los grupos (véase los ejercicios 4 y 7).
Pregunta 2: ¿Existe alguna relación o correlación entre las variables del estudio? Esta pregunta nos permitirá conocer la relación entre las variables de estudio para saber si puede existir cierta redundancia en la medida de la anatomía de la columna lumbar. Además, nos permitirá saber si existe un patrón significativo de influencia de una variable determinada a otra/otras variable/s.
Esta cuestión podrá resolverse mediante los estudios de regresión líneal y regresión múltiple (véase los ejercicios 5 y 6). Además, se realizará un análisis de clústering (véase ejercicio 7) para observar si es posible agrupar los datos de nuestra muestra.
El primer elemento de nuestros datos que analizaremos es la estadística básica y el comportamiento de cada variable entre los grupos.
Primeramente observaremos la estadística básica de cada variable mediante el comando “summary()”. Éste nos muestra los valores mínimos, máximos, primer y tercer cuartil, mediana y media para cada variable. Más adelante calcularemos su desviación estándar para cada grupo.
summary(ortho_df)
## pelvic_incidence pelvic_tilt lumbar_lordosis_angle sacral_slope
## Min. : 26.15 Min. :-6.555 Min. : 14.00 Min. : 13.37
## 1st Qu.: 46.43 1st Qu.:10.667 1st Qu.: 37.00 1st Qu.: 33.35
## Median : 58.69 Median :16.358 Median : 49.56 Median : 42.40
## Mean : 60.50 Mean :17.543 Mean : 51.93 Mean : 42.95
## 3rd Qu.: 72.88 3rd Qu.:22.120 3rd Qu.: 63.00 3rd Qu.: 52.70
## Max. :129.83 Max. :49.432 Max. :125.74 Max. :121.43
## pelvic_radius degree_spondylolisthesis class_2
## Min. : 70.08 Min. :-11.058 Normal :100
## 1st Qu.:110.71 1st Qu.: 1.604 Abnormal:210
## Median :118.27 Median : 11.768
## Mean :117.92 Mean : 26.297
## 3rd Qu.:125.47 3rd Qu.: 41.287
## Max. :163.07 Max. :418.543
## class_3
## Normal :100
## Hernia : 60
## Spondylolisthesis:150
##
##
##
Obsérvese que podemos determinar los límites anotados para cada variable y el comportamiento general mediante la media. Estos valores nos servirán más adelante para determinar un intervalo para presentar las variables en gráficos.
| Variable | Características |
|---|---|
| “pelvic_incidence” | valores de 26.15º hasta 129.83º |
| “pelvic_tilt” | valores de -6.555º hasta 49.432º |
| “lumbar_lordosis_angle” | valores de 14º hasta 125.74º |
| “sacral_slope” | valores de 13.37º hasta 121.43º |
| “pelvic_radius” | valores desde 70.08º hasta 163.07º |
| “degree_spondylolisthesis” | valores desde -11.058 mm hasta 418.543 mm |
| “class_2 y class_3” | incluyen 210 personas en el grupo “Abnormal” (60 con hernia y 150 con espondilolistesis) y 100 personas en el grupo “Normal”. |
Aunque esta descripción nos sirve para acotar los valores de las variables con los que trabajamos, los datos comprenden todos grupos de estudio en conjunto. Nuestro objetivo principal es la comparación entre personas del grupo “Normal” (sin anomalías aparentes relativas a la postura pélvica y lumbar) y “Abnormal” (personas con hernia o espondilolistesis diagnosticada) para determinar si existen diferencias significativas en las variables.
El hecho de dividir los datos podría llegar a mostrar una diferencia en el comportamiento de los mismos. Es por eso que realizaremos un análisis de comparación de variables según los grupos. El valor que compararemos será la media de cada variable, con su desviación estándar anotada. Asimismo, para hacer un análisis visual, las diferencias entre grupos de estudio se representarán en gráficos de cajas (boxplots). En cada caso se utilizarán los intervalos calculados anteriormente para que la representación de cada gráfico esté en el mismo intervalo y que así su comparación sea más intuitiva y fácil.
Hace falta remarcar que este análisis lo haríamos para las variables de tipo angular (las cinco primeras variables). No obstante, al tratarse de un trabajo de ámbito educativo, para la variable “grado de espondilolistesis” realizaremos la comparación igualmente, pero de antemano sabemos que será superior en pacientes del grupo “espondilolistesis” respecto al resto, ya que tienen dicha afectación diagnosticada.
Mediante la herramienta “aggregate” calculamos la media y la desviación estándar para cada grupo de nuestras variables.
aggregate(ortho_df[, 1:6], list(ortho_df$class_3), mean)
## Group.1 pelvic_incidence pelvic_tilt lumbar_lordosis_angle
## 1 Normal 51.68524 12.82141 43.54260
## 2 Hernia 47.63841 17.39879 35.46352
## 3 Spondylolisthesis 71.51422 20.74804 64.11011
## sacral_slope pelvic_radius degree_spondylolisthesis
## 1 38.86383 123.8908 2.186572
## 2 30.23961 116.4750 2.480251
## 3 50.76619 114.5188 51.896687
aggregate(ortho_df[, 1:6], list(ortho_df$class_3), sd)
## Group.1 pelvic_incidence pelvic_tilt lumbar_lordosis_angle
## 1 Normal 12.36816 6.778503 12.361388
## 2 Hernia 10.69713 7.016708 9.767795
## 3 Spondylolisthesis 15.10934 11.506169 16.397068
## sacral_slope pelvic_radius degree_spondylolisthesis
## 1 9.624004 9.014246 6.307483
## 2 7.555388 9.355720 5.531177
## 3 12.318813 15.579995 40.108030
Representaremos estos estadísticos a continuación para cada variable.
Variable 1: Incidencia pélvica (pelvic_incidence)
par(mfrow=c(1,4))
boxplot(ortho_df$pelvic_incidence[ortho_df$class_3=="Normal"],col="darkolivegreen1",main="Incidencia pélvica", xlab="Normal", ylim=c(20,130))
boxplot(ortho_df$pelvic_incidence[ortho_df$class_2=="Abnormal"],col="Tomato",main="Incidencia pélvica", xlab="Abnormal",ylim=c(20,130))
boxplot(ortho_df$pelvic_incidence[ortho_df$class_3=="Hernia"],col="coral1",main="Incidencia pélvica", xlab="Hernia", ylim=c(20,130))
boxplot(ortho_df$pelvic_incidence[ortho_df$class_3=="Spondylolisthesis"],col="darkturquoise",main="Incidencia pélvica", xlab="Spondylolisthesis", ylim=c(20,130))
par(mfrow=c(1,1))
Observamos que la media del grupo “Normal” es ligeramente inferior a la del grupo “Abnormal”. Al haber más pacientes en la muestra con “Spondylolisthesis” que en la de “Hernia”, la media del grupo “Abnormal” resulta ser superior a la del grupo “Normal”, con una desviación estándar muy alta. Este hecho se observa en los gráficos, donde parece que el grupo “Spondylolisthesis” supera significativamente la media otorgada por el grupo “Normal”. Como sólo observamos tres outliers aparentes en el grupo de espondilolistesis (la muestra no parece tener mucha distorsión), nuestra primera hipótesis será que los pacientes con espondilolistesis podrían tener una mayor incidencia pélvica general. Así, es probable que la variable “incidencia pélvica” tenga cierto peso a la hora de analizarse en dicho grupo.
Variable 2: Inclinación pélvica (pelvic_tilt)
par(mfrow=c(1,4))
boxplot(ortho_df$pelvic_tilt[ortho_df$class_3=="Normal"],col="darkolivegreen1",main="Inclinación pélvica", xlab="Normal", ylim=c(-7,50))
boxplot(ortho_df$pelvic_tilt[ortho_df$class_2=="Abnormal"],col="Tomato",main="Inclinación pélvica", xlab="Abnormal",ylim=c(-7,50))
boxplot(ortho_df$pelvic_tilt[ortho_df$class_3=="Hernia"],col="coral1",main="Inclinación pélvica", xlab="Hernia", ylim=c(-7,50))
boxplot(ortho_df$pelvic_tilt[ortho_df$class_3=="Spondylolisthesis"],col="darkturquoise",main="Inclinación pélvica", xlab="Spondylolisthesis", ylim=c(-7,50))
par(mfrow=c(1,1))
Se observa que la media del grupo “Normal” para la inclinación pélvica es ligeramente inferior que la de los otros grupos. No obstante, al existir una desviación estándar realmente alta para todos los grupos en esta variable (es decir, existe una gran variabilidad), es probable que la distinción entre grupos no llegue a ser significativa. Por lo tanto, no esperamos diferencias significativas entre individuos de los distintos grupos para esta variable.
Variable 3: Ángulo de lordosis lumbar (lumbar_lordosis_angle)
par(mfrow=c(1,4))
boxplot(ortho_df$lumbar_lordosis_angle[ortho_df$class_3=="Normal"],col="darkolivegreen1",main="Lordosis lumbar", xlab="Normal", ylim=c(13,126))
boxplot(ortho_df$lumbar_lordosis_angle[ortho_df$class_2=="Abnormal"],col="Tomato",main="Lordosis lumbar", xlab="Abnormal",ylim=c(13,126))
boxplot(ortho_df$lumbar_lordosis_angle[ortho_df$class_3=="Hernia"],col="coral1",main="Lordosis lumbar", xlab="Hernia", ylim=c(13,126))
boxplot(ortho_df$lumbar_lordosis_angle[ortho_df$class_3=="Spondylolisthesis"],col="darkturquoise",main="Lordosis lumbar", xlab="Spondylolisthesis", ylim=c(13,126))
par(mfrow=c(1,1))
En este caso observamos que la media del grupo “Normal” es parecida a la del grupo “Hernia”, pero difiere bastante de la del grupo “Spondylolisthesis”. Aunque exista algún outlier apartado que pueda hacer aumentar la desviación estándar de cada grupo (véase los puntos superiores que se distancian de la media en cada gráfico), es probable que esta comparativa pueda ser significativa. Así, nuestra segunda hipótesis es que los pacientes con espondilolistesis podrían tener un mayor ángulo de lordosis lumbar. Y por lo tanto, la variable “ángulo de lordosis lumbar” podría tener más peso en este grupo de pacientes.
Variable 4: Pendiente sacral (sacral_slope)
par(mfrow=c(1,4))
boxplot(ortho_df$sacral_slope[ortho_df$class_3=="Normal"],col="darkolivegreen1",main="Pendiente sacral", xlab="Normal", ylim=c(12,122))
boxplot(ortho_df$sacral_slope[ortho_df$class_2=="Abnormal"],col="Tomato",main="Pendiente sacral", xlab="Abnormal",ylim=c(12,122))
boxplot(ortho_df$sacral_slope[ortho_df$class_3=="Hernia"],col="coral1",main="Pendiente sacral", xlab="Hernia", ylim=c(12,122))
boxplot(ortho_df$sacral_slope[ortho_df$class_3=="Spondylolisthesis"],col="darkturquoise",main="Pendiente sacral", xlab="Spondylolisthesis", ylim=c(12,122))
par(mfrow=c(1,1))
Para la pendiente sacral observamos que las medias de los grupos difieren muy ligeramente, siendo superior en el grupo “Spondylolisthesis”. Existen outliers que pueden alterar la media de los datos, sobretodo en el grupo “Spondylolisthesis”, donde se observan varios puntos que se alejan significativamente del comportamiento del resto de las muestras. Es por eso que probablemente no existan diferencias significativas entre la media de los grupos para esta variable, ya que sin estos outliers la media del grupo “Spondylolisthesis” podría disminuir y se asemejaría aún más a la de los otros grupos.
Variable 5: Radio pélvico (pelvic_radius)
par(mfrow=c(1,4))
boxplot(ortho_df$pelvic_radius[ortho_df$class_3=="Normal"],col="darkolivegreen1",main="Radio pélvico", xlab="Normal", ylim=c(69,164))
boxplot(ortho_df$pelvic_radius[ortho_df$class_2=="Abnormal"],col="Tomato",main="Radio pélvico", xlab="Abnormal",ylim=c(69,164))
boxplot(ortho_df$pelvic_radius[ortho_df$class_3=="Hernia"],col="coral1",main="Radio pélvico", xlab="Hernia", ylim=c(69,164))
boxplot(ortho_df$pelvic_radius[ortho_df$class_3=="Spondylolisthesis"],col="darkturquoise",main="Radio pélvico", xlab="Spondylolisthesis", ylim=c(69,164))
par(mfrow=c(1,1))
Las medias de todos los grupos son semejantes, siendo ligeramente inferiores aquellas del grupo abnormal (incluyendo a los grupos hernia y espondilolistesis por separado). La gran variabilidad en el grupo de espondilolistesis (fíjese en su gran desviación estándar) hace dudar de la significación de esta comparación, aunque sea posible que la variable “radio pélvico” sea distintiva en el grupo “Normal” respecto al resto. Aún así, debido a la variabilidad no esperamos que se cumpla esta hipótesis a priori.
Variable 6: Grado de espondilolistesis (degree_spondylolisthesis)
par(mfrow=c(1,4))
boxplot(ortho_df$degree_spondylolisthesis[ortho_df$class_3=="Normal"],col="darkolivegreen1",main="G.espondilolistesis", xlab="Normal", ylim=c(-12,419))
boxplot(ortho_df$degree_spondylolisthesis[ortho_df$class_2=="Abnormal"],col="Tomato",main="G.espondilolistesis", xlab="Abnormal",ylim=c(-12,419))
boxplot(ortho_df$degree_spondylolisthesis[ortho_df$class_3=="Hernia"],col="coral1",main="G.espondilolistesis", xlab="Hernia", ylim=c(-12,419))
boxplot(ortho_df$degree_spondylolisthesis[ortho_df$class_3=="Spondylolisthesis"],col="darkturquoise",main="G.espondilolistesis", xlab="Spondylolisthesis", ylim=c(-12,419))
par(mfrow=c(1,1))
Efectivamente el grado de espondilolistesis en pacientes con espondilolistesis es superior respecto al resto, el cual varía ligeramente de la media 0. Esta comparación hace clara la presencia de un elemento de la muestra probablemente mal tomado (existe un paciente con 400 mm (40 cm) de desplazamiento vertebral, hecho bastante cuestionable). No obstante, la comparativa sigue vigente aun ignorar los outliers que distorsionan la media.
Por otro lado, fíjese que esta variable es la única que se mide en milímetros, mostrándonos que la media de variación de la distancia entre vértebras lumbares en pacientes con espondilolistesis es de 51.89 mm. Por lo tanto, la mayoría de pacientes con espondilolistesis del estudio tienen un grado de espondilolistesis positivo, es decir, la vértebra desplazada tiene un desplazamiento anterior respecto su vértebra inferior en la mayoría de casos.
Para realizar la regresión lineal eligiremos las dos variables que parezcan más relacionadas de entre todas las disponibles, siempre que esa relación analizada tenga un sentido real biológico. Para ver el comportamiento que cada variable tiene con las otras variables del estudio, utilizaremos un gráfico de dispersión de pares de variables (pairs). Nótese que eliminaremos del análisis las dos variables de tipo caracter (class_2 y class_3) debido a que no son numéricas. Para analizar la relación entre variables de forma cuantitativa, calcularemos las correlaciones entre las variables numéricas mediante el comando “cor()”.
pairs(ortho_df[, -c(7:8)])
cor(ortho_df[, -c(7:8)])
## pelvic_incidence pelvic_tilt lumbar_lordosis_angle
## pelvic_incidence 1.0000000 0.62919877 0.71728236
## pelvic_tilt 0.6291988 1.00000000 0.43276386
## lumbar_lordosis_angle 0.7172824 0.43276386 1.00000000
## sacral_slope 0.8149600 0.06234529 0.59838689
## pelvic_radius -0.2474672 0.03266781 -0.08034361
## degree_spondylolisthesis 0.6387427 0.39786228 0.53366701
## sacral_slope pelvic_radius degree_spondylolisthesis
## pelvic_incidence 0.81495999 -0.24746721 0.63874275
## pelvic_tilt 0.06234529 0.03266781 0.39786228
## lumbar_lordosis_angle 0.59838689 -0.08034361 0.53366701
## sacral_slope 1.00000000 -0.34212835 0.52355746
## pelvic_radius -0.34212835 1.00000000 -0.02606501
## degree_spondylolisthesis 0.52355746 -0.02606501 1.00000000
# Además, para que sea más facil interpretar las relaciones entre las variables vamos a crear una visualización de las interacciones en contexto de clasificación Normal-Anormal y Normal-Hernia-Espondilolistesis.
model.matrix(~0. + pelvic_incidence+pelvic_tilt+pelvic_radius+sacral_slope+lumbar_lordosis_angle+ degree_spondylolisthesis + class_2, ortho_df) %>%
cor() %>%
ggcorrplot(lab=T,
lab_size=2.5,
type = "lower")
model.matrix(~0. + pelvic_incidence+pelvic_tilt+pelvic_radius+sacral_slope+lumbar_lordosis_angle+ degree_spondylolisthesis + class_3, ortho_df) %>%
cor() %>%
ggcorrplot(lab=T,
lab_size=2.5,
type = "lower")
El comando “cor()” nos muestra los coeficientes de correlación entre cada variable, por defecto mediante el método de Pearson. Cuando más próximos a 1 o -1, más relación tienen las dos variables, y por lo tanto, el incremento (en caso de coeficiente positivo) o el decremento (en caso de coeficiente negativo) de una variable influencia al incremento/decremento de la otra variable, respectivamente.
Elección de variables De todos los coeficientes de correlación de Pearson entre dos variables distintas, el que tiene la cifra más elevada en valor absoluto es el de la incidencia pélvica (pelvic_incidence) con la pendiente sacral (sacral_slope), con un coeficiente de 0.8149600. Por lo tanto, para la regresión lineal, sería lógico analizar estas dos variables, que son las que tienen un mayor coeficiente de correlación. No obstante, según la bibliografía, la incidencia pélvica (PI) se calcula mediante la suma de la pendiente sacral (SS) más la inclinación pélvica (PT). Véase su descripción y definición indicada en el artículo de J.C.Le Huec et al. (2011):
https://www.ncbi.nlm.nih.gov/pmc/articles/PMC3175921/
Por lo tanto, no tendría mucho sentido analizar la relación entre estas dos variables, ya que son dos variables que se sabe que están relacionadas desde un principio. Y lo mismo pasaría con la relación entre la incidencia pélvica con la inclinación pélvica. Es por eso que decidiremos comprobar la relación entre dos variables diferentes, pero que conserven un coeficiente de correlación elevado, como el ángulo de lordosis lumbar (lumbar_lordosis_angle) con la incidencia pélvica (pelvic_incidence), que tiene un coeficiente de correlación distintivo de 0.7172824.
Así, en este ejercicio queremos comprobar si es cierto que el ángulo de lordosis lumbar es una variable que se correlaciona con la incidencia pélvica. Al tener un coeficiente de correlación relativamente elevado (es próximo a 1), nuestra hipótesis es que es probable que sí e cumpla esta relación y que ésta será positiva.
De las comparaciones visuales anteriores, si nos fijamos en el gráfico de comparación entre estas dos variables, observamos cómo la nube de puntos de los datos de nuestra muestra forma una línea en diagonal ascendente. Este hecho podría intuir que sí que existe una relación entre variables y verifica que éste tenga un valor positivo: cuando una de las variables aumenta, podría ser que la otra también. A continuación vamos a comprobar estadísticamente si este comportamiento es verídico o no.
Modelo de regresión lineal Para determinar si la relación es significativa realizaremos un análisis de regresión lineal entre estas dos variables. Para ello utilizaremos un modelo de regresión lineal (lm) que denominaremos como “ortho_lm”. Nosotros definiremos la incidencia pélvica como la variable independiente (explicativa) y el ángulo de lordosis lumbar como la variable dependiente (explicada), ya que el hecho de tener una incidencia pélvica más elevada podría implicar que el sacro de la persona en cuestión se inclinara hacia atrás y por lo tanto ésta tendría una modificación en su ángulo de lordosis lumbar.
ortho_lm <-lm(ortho_df$lumbar_lordosis_angle~ortho_df$pelvic_incidence)
summary(ortho_lm)
##
## Call:
## lm(formula = ortho_df$lumbar_lordosis_angle ~ ortho_df$pelvic_incidence)
##
## Residuals:
## Min 1Q Median 3Q Max
## -57.083 -7.991 -1.266 7.815 75.100
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 5.22081 2.68806 1.942 0.053 .
## ortho_df$pelvic_incidence 0.77211 0.04274 18.066 <2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 12.95 on 308 degrees of freedom
## Multiple R-squared: 0.5145, Adjusted R-squared: 0.5129
## F-statistic: 326.4 on 1 and 308 DF, p-value: < 2.2e-16
Mediante el “summary” del modelo observamos distintos valores estadísticos. El valor principal que nos interesa es el coeficiente de determinación (el R-squared) del modelo, que nos dictará qué tan cerca están los datos de una regresión lineal. En este caso, el coeficiente de determinación es relativamente próximo a 1 o -1 (siendo 0.5145) por lo que la bondad del modelo no parece demasiado fuerte, pero al ser superior a 0.5 podría corroborar nuestra hipótesis. En otras palabras, según este coeficiente, este modelo explicaría un 51.45% de la variabilidad de los datos.
Normalidad y varianza de los errores Para comprobar si el coeficiente obtenido es válido, analizaremos la normalidad de los errores del modelo mediante un histograma básico, un boxplot y el gráfico de cuartiles de los residuos.
residuos <-rstandard(ortho_lm)
par(mfrow=c(1,3))
hist(residuos, main = "Histograma de residuos")
boxplot(residuos, main = "Diagrama de cajas de residuos")
qqnorm(residuos, main = "G. cuartiles de residuos")
qqline(residuos)
par(mfrow=c(1,1))
En el primer gráfico observamos cómo, a nivel visual, los residuos parecen tener un comportamiento normal (la forma de una pendiente simétrica con un pico próximo al cero). El segundo gráfico nos añade la información relativa a los outliers. Parece ser que existen diversos valores de nuestra muestra que hacen alejarse los errores de la normalidad de los datos. No obstante, según el tercer gráfico, no parece ser significativo, ya que los residuos siguen una recta diagonal con una desviación en el tercer cuartil únicamente. Así, podemos concluir que los residuos parecen seguir una distribución normal.
Una vez verificada la normalidad de los residuos, analizaremos si la varianza de los errores es constante. El hecho de que la varianza de los errores sea constante implicaría que la variable dependiente variaría debido a nuestra variable independiente, no debido a otros factores externos a nuestro análisis. Para ello, compararemos el gráfico de nuestro modelo con los residuos estandarizados.
plot(fitted.values(ortho_lm),rstandard(ortho_lm),xlab="Valores ajustados",ylab="Residuos estandarizados", col="darkorchid3")
abline(h=0)
Con este gráfico verificamos que existen algunas anotaciones de nuestra muestra (outliers) que se alejan del comportamiento de la mayoría (es decir, que se alejan de la línea centrada al cero). Aun así, éstas parecen no ser significativas debido a que son pocas. A parte de estas anotaciones, el comportamiento de la nube de puntos del plot es constante y se adapta a la línea que marca el valor 0 de los residuos estandarizados.
Es por eso que aceptamos el modelo de regresión lineal y suponemos que el hecho de que su coeficiente de determinación no llegue a ser tan elevado como su coeficiente de correlación podría ser debido a la influencia de estos outliers.
Modelo final y conclusiones Finalmente, creemos que el modelo estudiado parece válido y que por lo tanto parece no existir una relación significativa entre las dos variables de estudio. Representaremos el modelo final en el siguiente gráfico.
plot(ortho_df$pelvic_incidence,ortho_df$lumbar_lordosis_angle, main="Modelo de regresión lineal", xlab="Incidencia pélvica", ylab= "Ángulo de lordosis lumbar")
abline(ortho_lm)
Como conclusión podemos decir que, como el coeficiente de determinación es relativamente elevado y el modelo de regresión es significativo, las dos variables analizadas (ángulo de lordosis lumbar e incidencia pélvica) parecen tener una relación en cuanto a su comportamiento. Esta relación es positiva (cuando el ángulo de una aumenta, el ángulo de la otra también). Este resultado es interesante, ya que ambas variables tienen el eje de medida en sitios distintos de la pelvis, a nivel anatómico. Aun así, los cambios en el comportamiento de una de las dos variables afectan a la rotación pélvica de forma que también se ocasiona cambios en el comportamiento de la otra variable.
Con el análisis del ejercicio 5 hemos podido responder parcialmente la pregunta 2 que nos proponíamos: ¿existe alguna relación entre nuestras variables de estudio? Mediante la regresión lineal hemos comprobado que existe una correlación positiva significativa entre el ángulo de lordosis lumbar y la incidencia pélvica. No obstante, también hemos visto unos índices de correlación elevados entre otras variables de estudio. Es por eso que ahora realizaremos una regresión múltiple para analizar qué variables tienen relación directa con otras.
Elección de variables Para definir un modelo inicial de regresión múltiple, debemos indagar de nuevo en los coeficientes de correlación de nuestras variables, para así elegir aquellos que sean próximos a 1 o -1. Nótese que trabajaremos con el mismo método de Pearson para calcularlos.
cor(ortho_df[, -c(7:8)])
## pelvic_incidence pelvic_tilt lumbar_lordosis_angle
## pelvic_incidence 1.0000000 0.62919877 0.71728236
## pelvic_tilt 0.6291988 1.00000000 0.43276386
## lumbar_lordosis_angle 0.7172824 0.43276386 1.00000000
## sacral_slope 0.8149600 0.06234529 0.59838689
## pelvic_radius -0.2474672 0.03266781 -0.08034361
## degree_spondylolisthesis 0.6387427 0.39786228 0.53366701
## sacral_slope pelvic_radius degree_spondylolisthesis
## pelvic_incidence 0.81495999 -0.24746721 0.63874275
## pelvic_tilt 0.06234529 0.03266781 0.39786228
## lumbar_lordosis_angle 0.59838689 -0.08034361 0.53366701
## sacral_slope 1.00000000 -0.34212835 0.52355746
## pelvic_radius -0.34212835 1.00000000 -0.02606501
## degree_spondylolisthesis 0.52355746 -0.02606501 1.00000000
En este caso observamos que, a parte de la correlación significativa entre la incidencia pélvica y la pendiente sacral, y entre la incidencia pélvica y el ángulo de lordosis lumbar, observados anteriormente, existen asociaciones con coeficientes de correlación elevados.
Es el caso de la incidencia pélvica con el grado de espondilolistesis (con un coeficiente de 0.6387427), seguido de la incidencia pélvica con la inclinación pélvica (coeficiente de 0.6291988). Finalmente existe un coeficiente de correlación de 0.59838689 entre el ángulo de lordosis lumbar y la pendiente sacral, un coeficiente de 0.53366701 entre el ángulo de lordosis lumbar y el grado de espondilolistesis, y un coeficiente de 0.52355746 entre la pendiente sacral y el grado de espondilolistesis.
A partir de estos casos, tenemos coeficientes de correlación en valor absoluto inferiores a 0.5 por lo que no los incluiremos a priori como destacables.
Por lo tanto, observamos que la variable “incidencia pélvica” es la que parece tener más relación con las otras, teniendo los coeficientes de correlación de Pearson más elevados para cuatro de las otras variables. Es por eso que parece intuitivo que esta variable sea nuestra variable objetivo para comparar con las otras. No obstante, como se ha mencionado en el ejercicio anterior, la incidencia pélvica es una variable de medida postural que se relaciona directamente con otras de nuestras variables de estudio y, por lo tanto, no sería muy interesante analizar su comportamiento en este ejercicio, ya que desde un principio se sabe que tendrá relación extremadamente alta con dos de las variables disponibles. Es por eso que eligiremos otra variable distinta que pueda presentarnos una comparación más interesante.
En este caso, hemos elegido el ángulo de lordosis lumbar. De momento, gracias a la regresión lineal, sabemos que tiene relación con la incidencia pélvica, pero queremos indagar sobre el comportamiento que tiene con las demás variables. Incluiremos todas las otras variables en el modelo, aunque sabemos que probablemente aquellas con coeficientes de correlación inferiores a 0.5 sean descartadas más adelante. En todo caso, realizaremos el modelo completo para que la elección se haga estadísticamente.
Modelo de regresión múltiple Realizamos el modelo de regresión del ángulo de lordosis lumbar con las cinco variables numéricas restantes.
ortho_modelomu<-lm(lumbar_lordosis_angle~pelvic_incidence+pelvic_tilt+sacral_slope+pelvic_radius+degree_spondylolisthesis, data=ortho_df)
summary(ortho_modelomu)
##
## Call:
## lm(formula = lumbar_lordosis_angle ~ pelvic_incidence + pelvic_tilt +
## sacral_slope + pelvic_radius + degree_spondylolisthesis,
## data = ortho_df)
##
## Residuals:
## Min 1Q Median 3Q Max
## -76.720 -7.415 -1.261 6.878 70.183
##
## Coefficients: (1 not defined because of singularities)
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) -10.55489 8.59928 -1.227 0.2206
## pelvic_incidence 0.76884 0.06996 10.989 <2e-16 ***
## pelvic_tilt -0.11410 0.09650 -1.182 0.2380
## sacral_slope NA NA NA NA
## pelvic_radius 0.14092 0.05911 2.384 0.0177 *
## degree_spondylolisthesis 0.05166 0.02555 2.022 0.0441 *
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 12.76 on 305 degrees of freedom
## Multiple R-squared: 0.5334, Adjusted R-squared: 0.5272
## F-statistic: 87.15 on 4 and 305 DF, p-value: < 2.2e-16
El “summary” del modelo nos informa de la significación de la correlación de nuestra variable objetivo (ángulo de lordosis lumbar) con cada otra variable. Obsérvese que los coeficientes significativos son únicamente con la incidencia pélvica (pelvic_incidence), con p-valor inferior a 2e-16, con el radio pélvico (pelvic_radius), con p-valor de 0.0177 y finalmente con el grado de espondilolistesis (degree_spondylolisthesis), con p-valor de 0.0441. Los asteriscos nos informan visualmente de que éstas son las únicas tres variables con correlaciones significativas para el ángulo de lordosis lumbar.
Además, nos otorga el coeficiente de determinación (R-squared) del modelo, siendo éste de 0.5334. Esto nos indica que la bondad del modelo de regresión múltiple es aceptable, ya que es superior a 0.5, aunque ésta no sea no muy elevada (sería mejor si estuviera más cerca de 1/-1). En otras palabras, un 53.34% de la variabilidad de la variable en cuestión (ángulo de lordosis lumbar) viene dada por las otras tres variables elegidas como variables con correlación significativa (incidencia pélvica, radio pélvico y grado de espondilolistesis).
Ahora que sabemos nuestras variables predictoras, perfeccionaremos el modelo de regresión, eliminando las otras variables que parecen no estar relacionadas directamente con nuestra variable objetivo.
ortho_modelomu2<-lm(lumbar_lordosis_angle~pelvic_incidence+pelvic_radius+degree_spondylolisthesis, data=ortho_df)
summary(ortho_modelomu2)
##
## Call:
## lm(formula = lumbar_lordosis_angle ~ pelvic_incidence + pelvic_radius +
## degree_spondylolisthesis, data = ortho_df)
##
## Residuals:
## Min 1Q Median 3Q Max
## -73.212 -7.786 -1.234 7.237 67.987
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) -7.62991 8.24113 -0.926 0.3553
## pelvic_incidence 0.72149 0.05740 12.569 <2e-16 ***
## pelvic_radius 0.12307 0.05719 2.152 0.0322 *
## degree_spondylolisthesis 0.05328 0.02553 2.087 0.0378 *
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 12.77 on 306 degrees of freedom
## Multiple R-squared: 0.5312, Adjusted R-squared: 0.5266
## F-statistic: 115.6 on 3 and 306 DF, p-value: < 2.2e-16
Confirmamos que el coeficiente de determinación del modelo es parecido al anterior (0.5312) y, como es relativamente alto, aceptamos este modelo final de regresión múltiple. En este sentido, si aceptamos el modelo, la variable “ángulo de lordosis lumbar” tendría gran parte de su variabilidad explicada por la incidencia pélvica, el radio pélvico y el grado de espondilolistesis. Como todos sus coeficientes de correlación son positivos, el modelo indicaría que un aumento de las tres variables independientes implicaría un aumento en el ángulo de lordosis lumbar. Como la significación de la correlación ente el ángulo de lordosis lumbar y la incidencia pélvica es mayor que en el resto de comparaciones, parece ser que gran parte de la variabilidad total de nuestra variable objetivo vendría dada por la incidencia pélvica.
A nivel anatómico estos datos tienen cierto sentido. El radio pélvico tiene su vértice de medida en el inicio del sacro. Si el ángulo de lordosis lumbar aumenta, la pelvis retrocede hacia atrás, haciendo que el radio aumente. Al moverse el sacro, pasaría algo parecido para la incidencia pélvica, aumentando su ángulo respecto el eje de la cabeza femoral. Y, finalmente, un aumento en ángulo de lordosis podría implicar un movimiento entre vértebras, definiendo la espondilolistesis en algunos de los pacientes de la muestra.
Para decidir si se puede utilizar el test ANOVA, tenemos que confirmar inicialmente si las variables siguen una distribución normal además de confirmar si no son homoscedásticas (es decir, tienen la misma varianza).
# Empezamos mirando las variables y calculando los valores de sus medias y errores estándar.
ortho_df2 %>% group_by(type_of_measure) %>% summarise("Mean of angle" = mean(angles), "Standart deviation" = sd(angles))
## # A tibble: 5 × 3
## type_of_measure `Mean of angle` `Standart deviation`
## <chr> <dbl> <dbl>
## 1 lumbar_lordosis_angle 51.9 18.6
## 2 pelvic_incidence 60.5 17.2
## 3 pelvic_radius 118. 13.3
## 4 pelvic_tilt 17.5 10.0
## 5 sacral_slope 43.0 13.4
ortho_df2 %>% summarise("Mean of displacement" = mean(degree_spondylolisthesis), "SD of displacement" = sd(degree_spondylolisthesis))
## # A tibble: 1 × 2
## `Mean of displacement` `SD of displacement`
## <dbl> <dbl>
## 1 26.3 37.5
# Creamos la visualización de la distribución de las variables con "denisity plot" y con "boxplot", marcando los valores de las medias. En la distribución normal el valor de la media y el valor de la mediana y moda coinciden. En cuanto a la visualización, parece que la variable "pelvic_radius" y "sacral_slope" podrían seguir una distribución normal. También miramos si el hecho de pertenecer a una otra clasificación puede afectar drásticamente al tipo de distribución. Se evidencia un cambio importante de comportamiento de la variable "degree_spondilolistesis" al comparar los grupos "Normal" y "Abnormal", además del grupo "Espondilolistesis" con el resto.
pt <- ortho_df2 %>% ggplot(aes(angles, fill = type_of_measure)) + geom_density(alpha = 0.2)
ptc2 <- ortho_df2 %>% ggplot(aes(angles, fill = type_of_measure)) + geom_density(alpha = 0.2)+ facet_grid(~class_2)
ptc3 <- ortho_df2 %>% ggplot(aes(angles, fill = type_of_measure)) + geom_density(alpha = 0.2)+ facet_grid(~class_3)
pp <- list(pt,ptc2,ptc3)
ggarrange(nrow = 3, ncol= 1, plotlist = pp)
pds <- ortho_df2 %>% ggplot(aes(degree_spondylolisthesis)) + geom_density(alpha = 0.2, fill = "darkgreen")
pds2 <- ortho_df2 %>% ggplot(aes(degree_spondylolisthesis, fill = class_2)) + geom_density(alpha = 0.2)+ facet_grid(~class_2)
pds3 <- ortho_df2 %>% ggplot(aes(degree_spondylolisthesis, fill = class_3)) + geom_density(alpha = 0.2)+ facet_grid(~class_3)
ppds <- list(pds,pds2,pds3)
ggarrange(nrow = 3, ncol= 1, plotlist = ppds)
# Evaluamos los valores de las medias y de los errores estándar según la clasificación:
# Segun class_2 (grupos "Normal" y "Abnormal")
ortho_df2 %>% group_by(type_of_measure, class_2) %>%
summarise("Mean of angle" = mean(angles), "Standart deviation" = sd(angles))
## `summarise()` has grouped output by 'type_of_measure'. You can override using the `.groups` argument.
## # A tibble: 10 × 4
## # Groups: type_of_measure [5]
## type_of_measure class_2 `Mean of angle` `Standart deviation`
## <chr> <fct> <dbl> <dbl>
## 1 lumbar_lordosis_angle Normal 43.5 12.4
## 2 lumbar_lordosis_angle Abnormal 55.9 19.7
## 3 pelvic_incidence Normal 51.7 12.4
## 4 pelvic_incidence Abnormal 64.7 17.7
## 5 pelvic_radius Normal 124. 9.01
## 6 pelvic_radius Abnormal 115. 14.1
## 7 pelvic_tilt Normal 12.8 6.78
## 8 pelvic_tilt Abnormal 19.8 10.5
## 9 sacral_slope Normal 38.9 9.62
## 10 sacral_slope Abnormal 44.9 14.5
ortho_df2 %>% group_by(class_2) %>%
summarise("Mean of displacement" = mean(degree_spondylolisthesis), "SD of displacement" = sd(degree_spondylolisthesis))
## # A tibble: 2 × 3
## class_2 `Mean of displacement` `SD of displacement`
## <fct> <dbl> <dbl>
## 1 Normal 2.19 6.28
## 2 Abnormal 37.8 40.6
# Según class_3 (grupos "Normal", "Hernia" y "Spondylolisthesis")
ortho_df2 %>% group_by(type_of_measure, class_3) %>%
summarise("Mean of angle" = mean(angles), "Standart deviation" = sd(angles))
## `summarise()` has grouped output by 'type_of_measure'. You can override using the `.groups` argument.
## # A tibble: 15 × 4
## # Groups: type_of_measure [5]
## type_of_measure class_3 `Mean of angle` `Standart deviation`
## <chr> <fct> <dbl> <dbl>
## 1 lumbar_lordosis_angle Normal 43.5 12.4
## 2 lumbar_lordosis_angle Hernia 35.5 9.77
## 3 lumbar_lordosis_angle Spondylolisthesis 64.1 16.4
## 4 pelvic_incidence Normal 51.7 12.4
## 5 pelvic_incidence Hernia 47.6 10.7
## 6 pelvic_incidence Spondylolisthesis 71.5 15.1
## 7 pelvic_radius Normal 124. 9.01
## 8 pelvic_radius Hernia 116. 9.36
## 9 pelvic_radius Spondylolisthesis 115. 15.6
## 10 pelvic_tilt Normal 12.8 6.78
## 11 pelvic_tilt Hernia 17.4 7.02
## 12 pelvic_tilt Spondylolisthesis 20.7 11.5
## 13 sacral_slope Normal 38.9 9.62
## 14 sacral_slope Hernia 30.2 7.56
## 15 sacral_slope Spondylolisthesis 50.8 12.3
ortho_df2 %>% group_by(class_3) %>%
summarise("Mean of displacement" = mean(degree_spondylolisthesis), "SD of displacement" = sd(degree_spondylolisthesis))
## # A tibble: 3 × 3
## class_3 `Mean of displacement` `SD of displacement`
## <fct> <dbl> <dbl>
## 1 Normal 2.19 6.28
## 2 Hernia 2.48 5.49
## 3 Spondylolisthesis 51.9 40.0
# Creamos la visualización de los valores de las medias y los errores estándar de las variables.
angles_ms <- ortho_df2 %>% ggplot(aes(y= angles, x = type_of_measure, color = type_of_measure)) +
stat_summary(fun = mean,
shape = 18,
size = 1)+
stat_summary(fun = mean,
fun.min = function(x) {mean(x) - sd(x)},
fun.max =function(x) {mean(x) + sd(x)},
geom = "errorbar")+
theme(axis.ticks.x = element_blank(),
axis.text.x = element_blank())
angles_ms2 <- ortho_df2 %>% ggplot(aes(y= angles, x = type_of_measure, color = type_of_measure)) +
stat_summary(fun = mean,
shape = 18,
size = 1)+
stat_summary(fun = mean,
fun.min = function(x) {mean(x) - sd(x)},
fun.max =function(x) {mean(x) + sd(x)},
geom = "errorbar")+
theme(axis.ticks.x = element_blank(),
axis.text.x = element_blank()) +
facet_wrap(~class_2)
angles_ms3 <- ortho_df2 %>% ggplot(aes(y= angles, x = type_of_measure, color = type_of_measure)) +
stat_summary(fun = mean,
shape = 18,
size = 1)+
stat_summary(fun = mean,
fun.min = function(x) {mean(x) - sd(x)},
fun.max =function(x) {mean(x) + sd(x)},
geom = "errorbar")+
theme(axis.ticks.x = element_blank(),
axis.text.x = element_blank()) +
facet_wrap(~class_3)
degree_ms <- ortho_df2 %>% ggplot(aes(y= degree_spondylolisthesis, x = 0)) +
stat_summary(fun = mean,
shape = 18,
size = 1,
color = "green")+
stat_summary(fun = mean,
fun.min = function(x) {mean(x) - sd(x)},
fun.max =function(x) {mean(x) + sd(x)},
color = "green",
geom = "errorbar")+
theme(axis.ticks.x = element_blank(),
axis.text.x = element_blank())
degree_ms2 <- ortho_df2 %>% ggplot(aes(y= degree_spondylolisthesis, x = class_2, color = class_2)) +
stat_summary(fun = mean,
shape = 18,
size = 1)+
stat_summary(fun = mean,
fun.min = function(x) {mean(x) - sd(x)},
fun.max =function(x) {mean(x) + sd(x)},
geom = "errorbar")+
theme(axis.ticks.x = element_blank(),
axis.text.x = element_blank()) +
facet_wrap(~class_2)
degree_ms3 <- ortho_df2 %>% ggplot(aes(y= degree_spondylolisthesis, x = class_3, color = class_3)) +
stat_summary(fun = mean,
shape = 18,
size = 1)+
stat_summary(fun = mean,
fun.min = function(x) {mean(x) - sd(x)},
fun.max =function(x) {mean(x) + sd(x)},
geom = "errorbar")+
theme(axis.ticks.x = element_blank(),
axis.text.x = element_blank()) +
facet_wrap(~class_3)
ggarrange(plotlist= list(angles_ms, angles_ms2, angles_ms3) , nrow= 3, ncol = 1, legend = "left")
## Warning: Removed 5 rows containing missing values (geom_segment).
## Warning: Removed 5 rows containing missing values (geom_segment).
## Warning: Removed 5 rows containing missing values (geom_segment).
## Warning: Removed 5 rows containing missing values (geom_segment).
## Warning: Removed 5 rows containing missing values (geom_segment).
## Warning: Removed 5 rows containing missing values (geom_segment).
ggarrange(plotlist= list(degree_ms, degree_ms2, degree_ms3) , nrow= 3, ncol = 1, legend = "left")
## Warning: Removed 1 rows containing missing values (geom_segment).
## Warning: Removed 1 rows containing missing values (geom_segment).
## Warning: Removed 1 rows containing missing values (geom_segment).
## Warning: Removed 1 rows containing missing values (geom_segment).
## Warning: Removed 1 rows containing missing values (geom_segment).
## Warning: Removed 1 rows containing missing values (geom_segment).
#Resumimos los datos en un gráfico de tipo "boxplot" para evaluar los cuartiles y su relación con los valores de las medias:
ortho_df2 %>% ggplot(aes(y= angles, x = type_of_measure)) +
geom_boxplot(alpha = 0.3, mapping = aes(fill = type_of_measure)) +
stat_summary(fun = mean,
shape = 18,
size = 1,
mapping = aes(color = type_of_measure)) +
theme(axis.ticks.x = element_blank(),
axis.text.x = element_blank())
## Warning: Removed 5 rows containing missing values (geom_segment).
ortho_df2 %>% ggplot(aes(y= angles, x = type_of_measure)) +
geom_boxplot(alpha = 0.3, aes(fill = type_of_measure)) +
stat_summary(fun = mean,
geom = "point",
shape = 18, size = 4,
aes(color = type_of_measure))+
facet_wrap(~class_2)+
theme(axis.ticks.x = element_blank(),
axis.text.x = element_blank())
ortho_df2 %>% ggplot(aes(y= angles, x = type_of_measure)) +
geom_boxplot(alpha = 0.3, aes(fill = type_of_measure)) +
stat_summary(fun = mean,
geom = "point", shape = 18,
size = 4,
aes(color = type_of_measure)) +
facet_wrap(~class_3)+
theme(axis.ticks.x = element_blank(),
axis.text.x = element_blank())
ortho_df2 %>% ggplot(aes(y= degree_spondylolisthesis,
x = 0)) +
geom_boxplot(fill = "red",
alpha = 0.3) +
stat_summary(fun = mean, shape = 18,
size =1)+
theme(axis.ticks.x = element_blank(),
axis.text.x = element_blank())
## Warning: Removed 1 rows containing missing values (geom_segment).
ortho_df2 %>% ggplot(aes(x =class_2,
y= degree_spondylolisthesis,
fill= class_2),
color = black) +
geom_boxplot(alpha = 0.3,
aes(fill = class_2)) +
stat_summary(fun = mean,
geom = "point",
shape = 18, size = 5,
aes(color = class_2))+
facet_wrap(~class_2)+
theme(axis.ticks.x = element_blank(),
axis.text.x = element_blank())
ortho_df2 %>% ggplot(aes(x = class_3, y= degree_spondylolisthesis), color = black) +
geom_boxplot(alpha = 0.3, aes(fill = class_3)) +
stat_summary(fun = mean,
geom = "point",
shape = 18,
size = 5,
aes(color = class_3)) +
facet_wrap(~class_3) +
theme(axis.ticks.x = element_blank(),
axis.text.x = element_blank())
# Normality tests: Ahora vamos a comprobar si la distribución de las variables realmente sigue una distribución normal (hipótesis 0). Para evaluarlo, como tenemos un abarcamiento muestral grande, vamos a utilizar el test de Kolmogorov-Smirnov, y su visualización con "QQ-plot".
# Kolmogorov-Smirnov
# Para realizar este test sobre el conjunto de data vamos a utilizar función map - una variante de lapply() dentro del paquete tidyverse (librería purr)
ortho_df %>%
select(where(is.numeric))%>%
colnames() %>%
set_names() %>%
map(~ lillie.test(ortho_df[,.])) %>%
map_dfr(., tidy, .id = "variable")
## # A tibble: 6 × 4
## variable statistic p.value method
## <chr> <dbl> <dbl> <chr>
## 1 pelvic_incidence 0.0708 7.33e- 4 Lilliefors (Kolmogorov-Smirnov) n…
## 2 pelvic_tilt 0.0829 2.37e- 5 Lilliefors (Kolmogorov-Smirnov) n…
## 3 lumbar_lordosis_angle 0.0687 1.25e- 3 Lilliefors (Kolmogorov-Smirnov) n…
## 4 sacral_slope 0.0578 1.43e- 2 Lilliefors (Kolmogorov-Smirnov) n…
## 5 pelvic_radius 0.0593 1.06e- 2 Lilliefors (Kolmogorov-Smirnov) n…
## 6 degree_spondylolisthesis 0.173 5.13e-25 Lilliefors (Kolmogorov-Smirnov) n…
# Vamos a mirar si el hecho de quitar el outlier (valor 418 del grado de espondilolistesis) nos afecta al tipo de distribución
ortho_df3 %>%
select(where(is.numeric))%>%
colnames() %>%
set_names() %>%
map(~ lillie.test(ortho_df[,.])) %>%
map_dfr(., tidy, .id = "variable") # Parece que no afecta, así que procedemos a seguir con el análisis con el data-frame completo.
## # A tibble: 6 × 4
## variable statistic p.value method
## <chr> <dbl> <dbl> <chr>
## 1 pelvic_incidence 0.0708 7.33e- 4 Lilliefors (Kolmogorov-Smirnov) n…
## 2 pelvic_tilt 0.0829 2.37e- 5 Lilliefors (Kolmogorov-Smirnov) n…
## 3 lumbar_lordosis_angle 0.0687 1.25e- 3 Lilliefors (Kolmogorov-Smirnov) n…
## 4 sacral_slope 0.0578 1.43e- 2 Lilliefors (Kolmogorov-Smirnov) n…
## 5 pelvic_radius 0.0593 1.06e- 2 Lilliefors (Kolmogorov-Smirnov) n…
## 6 degree_spondylolisthesis 0.173 5.13e-25 Lilliefors (Kolmogorov-Smirnov) n…
### QQplot
ggqqplot(ortho_df2, "angles", facet.by = "class_3", color = "type_of_measure")
ggqqplot(ortho_df2, "angles", facet.by = "class_2", color = "type_of_measure")
ggqqplot(ortho_df2, "degree_spondylolisthesis", color = "class_2", facet.by = "class_2")
ggqqplot(ortho_df2, "degree_spondylolisthesis", color = "class_3", facet.by = "class_3")
### Homogenidad - vamos a mirar si los datos son homogéneos
ortho_df %>%
select(where(is.numeric))%>%
colnames() %>%
set_names() %>%
map(~ fligner.test(ortho_df[,.] ~ class_2, data = ortho_df)) %>%
map_dfr(., tidy, .id = "variable")
## # A tibble: 6 × 5
## variable statistic p.value parameter method
## <chr> <dbl> <dbl> <dbl> <chr>
## 1 pelvic_incidence 18.9 1.40e- 5 1 Fligner-Killeen test of…
## 2 pelvic_tilt 14.0 1.80e- 4 1 Fligner-Killeen test of…
## 3 lumbar_lordosis_angle 25.3 4.98e- 7 1 Fligner-Killeen test of…
## 4 sacral_slope 16.4 5.11e- 5 1 Fligner-Killeen test of…
## 5 pelvic_radius 9.29 2.30e- 3 1 Fligner-Killeen test of…
## 6 degree_spondylolisthesis 89.8 2.58e-21 1 Fligner-Killeen test of…
ortho_df %>%
select(where(is.numeric))%>%
colnames() %>%
set_names() %>%
map(~ fligner.test(ortho_df[,.] ~ class_3, data = ortho_df)) %>%
map_dfr(., tidy, .id = "variable")
## # A tibble: 6 × 5
## variable statistic p.value parameter method
## <chr> <dbl> <dbl> <dbl> <chr>
## 1 pelvic_incidence 7.21 2.72e- 2 2 Fligner-Killeen test of…
## 2 pelvic_tilt 27.3 1.19e- 6 2 Fligner-Killeen test of…
## 3 lumbar_lordosis_angle 19.3 6.56e- 5 2 Fligner-Killeen test of…
## 4 sacral_slope 7.28 2.63e- 2 2 Fligner-Killeen test of…
## 5 pelvic_radius 22.9 1.07e- 5 2 Fligner-Killeen test of…
## 6 degree_spondylolisthesis 125. 6.04e-28 2 Fligner-Killeen test of…
Tras realizar varios tests de normalidad se confirma que la distribución de los datos NO es normal, así que procedemos a realizar tests ANOVA para conseguir resultados fiables con la idea en mente de que la distribución puede ser insuficientemente normal. El test de homogenidad de los datos sugiere que NO hay homogenidad en éstos.
# Realizoms el test de anova para diferentes variables.
anova_PI <- aov(pelvic_incidence~class_3, ortho_df)
anova_S <- aov(degree_spondylolisthesis~class_3, ortho_df)
anova_PR <- aov(pelvic_radius~class_3, ortho_df)
anova_LLA <- aov(lumbar_lordosis_angle~class_3, ortho_df)
# La relación de las medias entre grupos para la incidencia pélvica
par(mfrow=c(2,2))
plot(anova_PI)
par(mfrow=c(1,1))
plot(TukeyHSD(anova_PI))
# La relación de las medias entre grupos para el grado de espondilolistesis
par(mfrow=c(2,2))
plot(anova_S)
par(mfrow=c(1,1))
plot(TukeyHSD(anova_S))
# La relación de las medias entre grupos para ángulo del radio pélvico
par(mfrow=c(2,2))
plot(anova_PR)
par(mfrow=c(1,1))
plot(TukeyHSD(anova_PR))
# La relación de las medias entre gurpos para el ángulo de lordosis lumbar
par(mfrow=c(2,2))
plot(anova_LLA)
par(mfrow=c(1,1))
plot(TukeyHSD(anova_LLA))
Los tests anteriores sugieren que la diferencia entre los valores de las medias de algunas variables difiere de forma significativa entre algunos grupos (ver conclusiones en el ejercicio 8).
#Creamos un scatter plot para evaluar la presencia de algún tipo de aglomeración.
ortho_df %>% ggplot(aes(lumbar_lordosis_angle, pelvic_incidence, color = class_3)) +
geom_point()
ortho_df %>% ggplot(aes(lumbar_lordosis_angle, pelvic_incidence, color = class_2)) +
geom_point()
# Miramos si se puede identificar algún número de clúster
fviz_nbclust(ortho_df[,c(1,3)], kmeans, method = "wss")
fviz_nbclust(ortho_df[,c(1,3)], kmeans, method = "silhouette")
# Probamos de clasificar las valores segun la distancia entre los puntos del scatter plot para 2 y 3 clústeres
cluster2 <- kmeans(ortho_df[, c(1,3)], 2, nstart = 25)
table(cluster2$cluster, ortho_df$class_2)
##
## Normal Abnormal
## 1 16 116
## 2 84 94
summary(cluster2)
## Length Class Mode
## cluster 310 -none- numeric
## centers 4 -none- numeric
## totss 1 -none- numeric
## withinss 2 -none- numeric
## tot.withinss 1 -none- numeric
## betweenss 1 -none- numeric
## size 2 -none- numeric
## iter 1 -none- numeric
## ifault 1 -none- numeric
fviz_cluster(cluster2, data = ortho_df[, c(1,3)])
cluster2 # Se evidencia que la clasificación entre 2 clústeres es bastente complicada y la búsqueda de algomeraciones por distancia no ha sido muy eficaz. Se acierta solo en 61.4% de los casos.
## K-means clustering with 2 clusters of sizes 132, 178
##
## Cluster means:
## pelvic_incidence lumbar_lordosis_angle
## 1 76.19685 68.75039
## 2 48.85381 39.45807
##
## Clustering vector:
## [1] 1 1 1 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2
## [38] 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2
## [75] 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2
## [112] 2 2 2 2 2 2 2 2 1 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 1 1 2 2 2 2
## [149] 2 2 2 2 1 2 2 1 2 2 2 1 2 2 2 1 1 1 1 2 2 2 2 2 1 1 2 2 1 1 1 2 2 1 1 1 1
## [186] 2 1 2 1 2 1 1 2 1 2 1 1 2 2 2 2 1 2 1 1 1 1 1 1 2 1 1 1 1 2 1 1 1 1 1 1 1
## [223] 1 1 1 1 1 2 1 1 1 1 1 1 1 1 1 1 1 1 2 1 1 2 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1
## [260] 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1
## [297] 1 1 1 1 1 1 1 1 1 1 1 1 1 1
##
## Within cluster sum of squares by cluster:
## [1] 43253.08 33224.06
## (between_SS / total_SS = 61.4 %)
##
## Available components:
##
## [1] "cluster" "centers" "totss" "withinss" "tot.withinss"
## [6] "betweenss" "size" "iter" "ifault"
cluster3 <- kmeans(ortho_df[, c(1,3)], 3, nstart = 25)
table(cluster3$cluster, ortho_df$class_3)
##
## Normal Hernia Spondylolisthesis
## 1 38 12 67
## 2 58 48 13
## 3 4 0 70
summary(cluster3)
## Length Class Mode
## cluster 310 -none- numeric
## centers 6 -none- numeric
## totss 1 -none- numeric
## withinss 3 -none- numeric
## tot.withinss 1 -none- numeric
## betweenss 1 -none- numeric
## size 3 -none- numeric
## iter 1 -none- numeric
## ifault 1 -none- numeric
fviz_cluster(cluster3, data = ortho_df[, c(1,3)])
cluster3 # La clusterización en 3 clústeres se acierta en 74.1% de casos.
## K-means clustering with 3 clusters of sizes 117, 119, 74
##
## Cluster means:
## pelvic_incidence lumbar_lordosis_angle
## 1 63.58681 53.46109
## 2 43.99372 35.12170
## 3 82.14937 76.54268
##
## Clustering vector:
## [1] 3 3 3 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2
## [38] 2 2 2 2 2 2 2 2 2 2 1 2 2 1 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2
## [75] 2 2 2 2 2 2 2 2 2 2 2 2 2 1 2 1 1 2 2 2 2 2 1 1 2 1 2 2 1 2 2 2 2 2 2 2 2
## [112] 1 2 2 2 1 2 2 2 1 2 1 2 2 2 2 2 2 2 2 2 1 2 2 1 2 1 1 1 2 1 1 1 1 1 1 1 1
## [149] 1 1 1 1 1 1 1 3 1 1 1 3 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1
## [186] 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 3 1 1 1 1 1 3 1 1 1 1 1 1 1 3 3 1 1 1 3 3
## [223] 1 1 1 1 3 1 3 3 1 3 1 1 3 1 1 1 3 3 1 3 1 1 3 1 1 1 3 3 3 3 1 3 3 1 3 3 3
## [260] 3 3 1 1 3 3 1 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 1 3 3 1
## [297] 3 3 1 3 3 3 3 3 3 3 3 3 3 3
##
## Within cluster sum of squares by cluster:
## [1] 17597.78 12236.25 21400.69
## (between_SS / total_SS = 74.1 %)
##
## Available components:
##
## [1] "cluster" "centers" "totss" "withinss" "tot.withinss"
## [6] "betweenss" "size" "iter" "ifault"
# Otra forma de realizar el análisis de clúster en estos casos es la clusterización jerárquica:
cluster <- hclust(dist(ortho_df), method="complete")
## Warning in dist(ortho_df): NAs introduced by coercion
plot(cluster)
pltree(agnes(ortho_df, method = "complete"))
Por lo tanto, parece ser que la bondad de la realización de clústeres es superior al utilizar 3 clústeres con nuestros datos, respecto 2.
Mediante este estudio hemos podido resolver las dos preguntas principales hechas en el principio del trabajo.
Pregunta 1: ¿Existen diferencias significativas entre las variables angulares del grupo “Normal” respecto al grupo “Abnormal”?.
Según las medias aritméticas de nuetra muestra, parece existir una diferencia significativa entre algunas de las variables según los distintos grupos (Normal/Hernia/Espondilolistesis). Primeramente hemos visto cómo las medias de los pacientes con espondilolistesis para la incidencia pélvica y el ángulo de lordosis lumbar son superiores que en sujetos del grupo “Normal”. Este hecho denota que el desliz de las vértebras de los pacientes con espondilolistesis afecta a estas dos variables métricas más que a los otros sujetos. Por lo tanto, los resultados sugieren que la incidencia pélvica y el ángulo de lordosis lumbar podrían llegar a ser dos variables definitorias para la espondilolistesis. Para los pacientes con “Hernia”, no obstante, no había diferencias significativas respecto a los sujetos del grupo “Normal”. Por lo tanto, el desplazamiento vertebral (que crea la espondilolistesis) influye claramente más que la hernia en las variables métricas posturales.
Por otro lado, tenemos que no parecen existir diferencias significativas entre los grupos para el radio pélvico, la pendiente sacral y la inclinación pélvica, por lo que éstas podrían ser variables menos explicativas para las afectaciones estudiadas. De hecho, el grupo de pacientes con Hernia no difiere significativamente del grupo Normal para ninguna variable, dando a entender que las variables métricas estudiadas en este análisis no serían buenos distintivos para esta patología. Quizás con más estudios con muestras mayores y más representativas se podría contrastar los datos de las dos patologías y verificar si las variables utilizadas en estos ejercicios son buenas variables tanto definidoras como métricas.
Finalmente, aunque creemos estas relaciones son lógicas y coherentes para la muestra, la fiabilidad del test de ANOVA es dudosa debido al tipo de distribución de sus datos (no parecen seguir una distribución normal). Por lo tanto no podemos verificar estadísticamente que se cumplan dichas relaciones con este test, y quizás sería preferible utilizar otros tests de comparación para determinarlo. Es remarcable que existen ciertos outliers en los datos de la muestra que podrían afectar al cómputo final de las comparaciones. Uno de estos outliers comprende los valores de un paciente con ángulos y desplazamientos vertebrales ciertamente anómalos para distintas variables, sugeriendo que es un error puntual de toma de muestras. En todo caso, al tratarse de un solo paciente, no creemos que sea decisivo para el comportamiento general de los datos en conjunto.
A modo de resumen, las variables distintivas entre grupos son:
| Variable | Relación |
|---|---|
| pelvic_incidence | NO varía entre el grupo Normal y el grupo Hernia, pero varía en comparación con el grupo Espondilolistesis |
| degree_spondylolisthesis | NO varía entre el grupo Normal y el grupo Hernia, pero varía en comparación con el grupo Espondilolistesis |
| pelvic_radius | NO varía entre el grupo Espondilolistesis y el grupo Hernia, pero varía en comparación con el grupo espondilolistesis y el grupo Normal |
| lumbar_lordosis_angle | Varía de forma importante entre los tres grupos. |
Pregunta 2: ¿Existe alguna relación o correlación entre las variables del estudio?
Según los resultados de nuestros análisis, parece existir una relación entre algunas de las variables.
Primeramente, mediante los coeficientes de correlación de Pearson, hemos detectado una posible correlación significativa entre la incidencia pélvica con la pendiente sacral y con la inclinación pélvica. Estas correlaciones no se han verificado mediante ningún análisis estadístico adicional debido a que han sido documentadas en la bibliografía anteriormente (J.C.Le Huec et al. (2011)).
En referencia al ángulo de lordosis lumbar, mediante un análisis de regresión múltiple hemos detectado su relación con incidencia pélvica, el radio pélvico y el grado de espondilolistesis. Adicionalmente, gracias al análisis de regresión lineal hemos podido observar más detalladamente la relación de este ángulo de lordosis lumbar con la incidencia pélvica y hemos comprobado que ésta parece ser más significativa para la incidencia pélvica que para las otras dos variables. Sin embargo, en todos los tres casos es una relación positiva. Es decir, según el modelo de estudio, un aumento de alguna de las tres variables predictoras (incidencia pélvica, radio pélvico o grado de espondilolistesis) podría provocar un aumento en el ángulo de lordosis lumbar.
Estos resultados tienen sentido a nivel anatómico y nos permiten determinar que el ángulo de lordosis lumbar es una variable importante a tener en cuenta en los estudios de medida pélvica.
Finalmente, durante el análisis de clúster se han realizado dos tipos de análisis (uno no jerárquico, utilizando valores de distancia entre puntos, y otro jerárquico). Al realizar el test “kmeans()” y al evaluar el data-frame se evidencia que el número óptimo de clústeres se sitúa entre 2 y 3, por lo que se ha realizado un análisis para las dos situaciones, objetivando mejor la calidad de aciertamiento de clusterización (un 74%) para la situación con 3 clústeres.